function sim = solve_bf(sim, eq_SS, param, glob, options)

% Solves "backwards" and iterates "forwards" on one iteration of the MIT
% shock

% start at the end and solve backwards
v2    = reshape(eq_SS.v.v2, glob.n(1), glob.n(2)); % assume that the expected value is in SS at date T

tic() 

for t = options.T:-1:1
    
    param.C = sim.C(t);
    param.A = sim.A(t); 
    param.w = sim.W(t);
    
    
    param.mu_f = sim.cF(t);
    param.ce = sim.cE(t);
    
 %   v2    = griddedInterpolant(glob.B, glob.S,v2, 'linear');   
    
 param.sdf = sim.sdf(t);
 
 
    %% find D that sets value of entry equal to 0
    tosolve = @(D) solve_b_minimal(D, v2, param, glob, options);
   
    if t == 1
        disp('t = 1')
    end
    
    %D = fzero(tosolve, eq_SS.D);
    try
        D = bisect(tosolve, .5*eq_SS.D, 2*eq_SS.D);
    catch
        disp('D not bracketed')
    end
        
    param.D = D;
    
    [~,ve, v] = solve_b(D, v2, param, glob, options);
    sim.ve(t) = ve.eve;
    sim.D(t) = D;
    % update v2
    v2    = reshape(v.v2, glob.n(1), glob.n(2));

    %% continue with incumbent problem
%     % solve on finer grid
%     tmp      = glob.PnK;
%     glob.PnK = glob.PnKf;
%     vf       = solve_valfunc([v2; v2],glob.sf,param,glob,options);
%     glob.PnK = tmp;
    
    vf = v;
    
    % Record optimal policy
    sim.y(:,t)  = vf.y; 
    sim.l(:,t)  = vf.l;
    sim.p(:,t)  = vf.p;  
    
    sim.exit(:,t) = vf.p_exit;
        
end
toc()

%% Now iterate forwards on policies
tic()

fspaceerg   = fundef({'spli',glob.bgridf,0,1});


for t = 1:options.T
    
    % Distribution from last period
   if t > 1
       L               = sim.mu(:,t-1);
       
       lp = sim.l(:,t-1);
       le = glob.ve_l;
       
       p_exit          = sim.exit(:,t-1);
   else
        L               = eq_SS.L;
        
        lp = eq_SS.l;
        le = eq_SS.ve_l;
        
        p_exit          = eq_SS.p_exit;
   end
       
    % distribution of entrants
    G                  = zeros(glob.Nsf,1);
    G((glob.sf(:,1) == glob.bgridf(end))) = glob.a_stat_E;

    entrant_policy = repmat(le, 1, glob.nf(1))';
    entrant_policy = entrant_policy(:)';
    
    entrant_policy          = min(entrant_policy, glob.maxb);
    entrant_policy          = max(entrant_policy, glob.minb);

    Qb             = funbas(fspaceerg, entrant_policy');

    Qent           = dprod(glob.QeI,Qb);
    G = Qent'*G;
   
    % set up transition matrix
    lp          = min(lp, glob.maxb);
    lp          = max(lp, glob.minb);
      
    Qb              = funbas(fspaceerg, lp);
    Q               = dprod(glob.Qef,Qb);
        
     param.D = sim.D(t);
     param.C = sim.C(t);
    param.w = sim.W(t);
    yf = sim.y(:,t);
    
    tosolve  = @(M) iterate_dist(M, L, G, Q, p_exit,yf, param,glob,options);
    
    try
        M        = bisect(tosolve, max(0,eq_SS.M - 1), eq_SS.M + 5);
    catch
        M = 0; 
    end
        
    sim.M(t) = M;
    
    
    sim.mu(:,t) = Q'*((1-param.gamma).*(1-p_exit).*L) + M*G;

    sim.entry_mass(t) = M;
    sim.entry_rate(t) = sim.entry_mass(t)/sum(sim.mu(:,t));

    eq_mit.L          = sim.mu(:,t);
    eq_mit.yf         = sim.y(:,t);
    
    aggs              = find_agg_mit(eq_mit, param, glob, options);
    
    sim.C(t)          = aggs.C;        
                                  
    L                 = sum(eq_mit.L.*sim.l(:,t));                    
    sim.L(t)          = L;
    sim.W(t)          = L^param.nu*param.psi;
      
    

    
end
toc()

for t = 1:options.T
        
    risk_aversion = param.ra;
    
    if t < options.T
        %param.sdf = param.beta*((sim.C(t+1) - param.psi/(1+param.nu)*sim.L(t+1)^(1+param.nu))^(-risk_aversion))/...
        %    ((sim.C(t) - param.psi/(1+param.nu)*sim.L(t)^(1+param.nu))^(-risk_aversion));
        param.sdf  = param.beta*(sim.C(t+1)/sim.C(t));
    else
        param.sdf = param.beta;
    end

    sim.sdf(t) = param.sdf;
    
end

end